เตรียมตัวก่อนเรียน
บทเรียนนี้จะกล่าวถึงการวิเคราะห์โมเดลพหุระดับด้วยวิธีการแบบเบส์ ซึ่งมีความแตกต่างไปจากการวิเคราะห์ด้วยวิธีการแบบดั้งเดิมที่เรียกว่าวิธีการแบบ frequentist ส่วนแรกของบทเรียนจะปูพื้นฐานเกี่ยวกับสถิติแบบเบส์ก่อน จากนั้นจึงกล่าวถึงการใช้สถิติแบบเบส์ในการวิเคราะห์โมเดลพหุระดับ
โปรแกรมที่ใช้ในการบรรยายประกอบด้วย
Jamovi
R
JAGS (Just Another Gibb Sampler)
STAN
ทั้งนี้เพื่อให้ไม่เป็นการเสียเวลา ขอให้ผู้เรียนดาวน์โหลดและติดตั้งโปรแกรมดังกล่าวให้พร้อมก่อนวันที่เรียน
ความรู้พื้นฐาน
หัวข้อนี้เป็นความรู้พื้นฐานเกี่ยวกับสถิติแบบดั้งเดิม และ สถิติแบบเบส์ โดยจะกล่าวถึงข้อตกลงเบื้องต้นของสถิติแต่ละประเภท และข้อจำกัดของสถิติแบบดั้งเดิมที่สามารถใช้สถิติแบบเบส์เข้ามาช่วยแก้ปัญหาได้ รายละเอียดมีดังนี้
สถิติแบบ Frequentist
สถิติแบบดั้งเดิมหรือที่เรียกกันว่า frequentist เป็นสถิติวิเคราะห์ที่ถูกใช้งานเป็นหลักในปัจจุบัน สถิติเกือบทุกตัวที่ผู้เรียนได้เรียนมาจนถึงปัจจุบันนี้ล้วนเป็นสถิติวิเคราะห์ภายใต้กระบวนทัศน์แบบ frequentist ทั้งสิ้น สถิติแบบ frequentist นี้มีข้อตกลงเบื้องต้นว่า (1) พารามิเตอร์เป็นค่าคงที่ ที่มีอยู่จริงในประชากรเป้าหมาย อย่างไรก็ตามด้วยข้อจำกัดที่ไม่สามารถเก็บรวบรวมข้อมูลจากทั้งประชากรได้ทำให้ผู้วิเคราะห์ไม่ทราบค่าพารามิเตอร์ดังกล่าว และจะประมาณโดยใช้ข้อมูลจากตัวอย่าง (2) ข้อมูลตัวอย่าง (sample data) เขียนแทนด้วย \(\bf{x}= \left \{\bf{x_1}, \bf{x_2},...,\bf{x_n} \right \}\) เป็นค่าที่สุ่มมาจากการประชากร ดังนั้นตัวอย่างที่สุ่มได้จึงเป็นค่าสุ่มที่และมีความเป็นไปได้ทั้งหมดเท่ากับ เมื่อ \(N\) คือขนาดประชากร และ \(n\) คือขนาดตัวอย่าง จากข้อตกลงเบื้องต้นข้อที่ (2) นี้จึงทำให้ (3) ค่าสถิติที่ประมาณได้จากข้อมูลตัวอย่างเป็นตัวแปรสุ่มที่มีการแจกแจงความน่าจะเป็น เรียกการแจกแจงความน่าจะเป็นของค่าสถิตินี้ว่า การแจกแจงความน่าจะเป็นของฟังก์ชันที่ได้จากตัวอย่างสุ่ม (sampling distribution)
ลักษณะการแจกแจงความน่าจะเป็นของค่าสถิติดังกล่าวมักประมาณได้ด้วย 2 วิธีการ วิธีการแรก ประมาณโดยใช้ทฤษฎีลิมิตลู่เข้าสู่ส่วนกลาง (central limit theorem) ซึ่งเป็นการคาดการณ์ลักษณะการแจกแจงของค่าสถิติในเชิงทฤษฎี เช่น การประมาณช่วงความเชื่อมั่นของค่าเฉลี่ยมีสูตรเป็น \(\overline{X}\pm t_{\alpha/2,n-1} \frac{S}{\sqrt{n}}\) จะเห็นว่าการประมาณช่วงความเชื่อมั่นดังกล่าวมีการใช้การแจกแจงที (student t distribution) ในการคำนวณส่วน margin of error การแจกแจงดังกล่าวพิสูจน์ในทางทฤษฎีโดยใช้ central limit theorem ได้ว่า เป็น sampling distribution ของ \(\overline{X}\) หากข้อมูลตัวอย่างที่ใช้ในการวิเคราะห์สุ่มมาจากประชากรที่มีการแจกแจงแบบปกติหรือใกล้เคียง แต่ไม่ทราบค่าความแปรปรวนของประชากร หรือถ้าหากประชากรไม่ได้มีการแจกแจงแบบปกติ แต่ตัวอย่างมีขนาดใหญ่เพียงพอ ก็ยังสามารถใช้ sampling distribution ดังกล่าวประมาณการแจกแจงของ \(\overline{X}\) ได้อยู่ หรือในการทดสอบสมมุติฐาน \(H_0: \mu=\mu_0\) สถิติทดสอบที่ใช้คือ \(t^*=\frac{\overline{X}-\mu_0}{S/ \sqrt{n}}\) จะมี sampling distributionในทางทฤษฎีเป็นการแจกแจงแบบที ที่มีองศาความเป็นอิสระเท่ากับ \(n-1\) เช่นเดียวกันหากข้อกำหนดเบื้องต้นของ central limit theorem เป็นจริงหรือใกล้เคียง การทราบการแจกแจงของค่าสถิติดังกล่าวทำให้ผู้วิเคราะห์สามารถประเมินความคลาดเคลื่อนในการตัดสินใจของการทดสอบสมมุติฐาน และสรุปผลได้ดังที่ได้เรียนไปแล้วในรายวิชาพื้นฐาน
ความถูกต้องของทฤษฎีนี้ขึ้นอยู่กับปัจจัยแวดล้อมที่เกี่ยวข้องว่ามีความสอดคล้องกับทฤษฎีมากน้อยเพียงใด ปัจจัยหนึ่งที่มีความสำคัญคือขนาดตัวอย่าง กล่าวคือหากขนาดตัวอย่างมีค่าน้อยเกินไป หรือการแจกแจงของข้อมูลมีความเบี่ยงเบนออกไปจากข้อตกลงเบื้องต้นของการวิเคราะห์ไปมาก ค่าสถิติที่คำนวณได้อาจไม่มีลักษณะการแจกแจงที่เป็นไปตามทฤษฎี และส่งผลให้การอนุมานเชิงสถิติมีความคลาดเคลื่อน อีกวิธีการหนึ่งที่สามารถใช้ในการประมาณ sampling distribution ของค่าสถิติได้คือใช้วิธีการสุ่มซ้ำ (resampling method) เช่น boostraping หรือ jackknife เพื่อประมาณแนวโน้มการแจกแจงของค่าสถิติ เป็นต้น
t.test(dat$Ach, mu=5.5)
##
## One Sample t-test
##
## data: dat$Ach
## t = 5.078, df = 149, p-value = 1.123e-06
## alternative hypothesis: true mean is not equal to 5.5
## 95 percent confidence interval:
## 5.709732 5.976934
## sample estimates:
## mean of x
## 5.843333
คำถาม
output ข้างต้นแสดงผลการประมาณช่วงความเชื่อมั่น และทดสอบสมมุติฐานทางสถิติเกี่ยวกับผลสัมฤทธิ์ทางการเรียนของนักเรียน (คะแนนเต็ม 10) โดยจากการประมาณด้วยช่วงความเชื่อมั่น 95% ในข้างต้นพบว่ามีค่าเท่ากับ \([5.71, 5.98]\) ซึ่งหมายความว่าอะไร?
จากผลการทดสอบสมมุติฐาน \(H_0: \mu_{Ach}=5.5\) พบว่ามีค่า p-value < .000 หมายความว่าอะไร?
เราสามารถตัดสินใจยอมรับสมมุติฐาน \(H_0\) ได้หรือไม่?
จากข้อตกลงเบื้องต้นของสถิติแบบ frequentist นี้ทำให้การวิเคราะห์ข้อมูลเกิดข้อจำกัดต่าง ๆ ดังนี้
ด้านที่สำคัญคือด้านการแปลผลการวิเคราะห์ที่จะเห็นว่าความหมายของผลการวิเคราะห์ต่าง ๆ เป็นการกล่าวถึงพารามิเตอร์ที่สนใจแบบอ้อม ๆ ทั้งสิ้น ไม่มีการวิเคราะห์ใดที่สามารถอ้างอิงไปยังพารามิเตอร์ที่สนใจได้โดยตรง ทั้งนี้เป็นเพราะการอนุมานเชิงสถิติแบบ frequentist พึ่งพาเครื่องมือหลักคือ samping distribution ของค่าสถิติ ในขณะที่พารามิเตอร์ถูกสมมุติให้เป็นค่าคงที่ตั้งแต่ต้น probability statement ต่าง ๆ จึงเป็นของค่าสถิติทั้งสิ้น ดังความความเชื่อมั่นที่รับประกันในช่วงความเชือมั่น และค่าความคลาดเคลื่อนในการทดสอบสมมุติฐานจึงเป็นค่าความน่าจะเป็นที่ใช้รับประกันความเป็นไปได้ของค่าสถิติทั้งสิ้นไม่ใช่ของค่าพารามิเตอร์ที่สนใจจริง ๆ ในการวิเคราะห์
ด้านที่สองคือด้านการประมาณค่าพารามิเตอร์ภายใต้โมเดลซับซ้อน เครื่องมือสำคัญสำหรับประมาณค่าพารามิเตอร์ในโมเดลต่าง ๆ ของสถิติแบบ frequentist คือตัวประมาณ เช่น ตัวประมาณกำลังสองน้อยสุด ตัวประมาณภาวะความควรจะเป็นสูงสุด ซึ่งล้วนเป็นวิธีการที่จะหาค่าประมาณที่ดีที่สุดโดยอิงจากจุดต่ำสุด หรือจุดสูงสุดของฟังก์ชันวัตถุประสงค์ ดังตัวอย่างในรูปด้านล่าง ซึ่งในกรณีที่โมเดลมีความซับซ้อนมาก ๆ ฟังก์ชันวัตถุประสงค์ดังกล่าวจะมีความซับซ้อนโดยมีจำนวนมิติที่สูงขึ้น และอาจมีจุดสูงต่ำ และ/หรือ จุดต่ำสุดสัมพัทธ์จำนวนมาก ซึ่งเป็นอุปสรรคให้การประมาณค่าด้วยวิธีการในข้างต้นอาจได้ค่าประมาณที่คลาดเคลื่อนไปจากค่าที่เหมาะสมที่สุด หรืออัลกอริทึมไม่สามารถหาค่าประมาณที่เหมาะสมได้ภายใต้จำนวนรอบของการประมาณที่กำหนด
- สืบเนื่องจากที่การอนุมานเชิงสถิติแบบ frequentist นั้นมักอิงกับทฤษฎีความน่าจะเป็นที่อิงข้อกำหนดเบื้องต้นที่ขนาดตัวอย่างต้องมีขนาดใหญ่เพียงพอ ซึ่งทำให้สถิติแบบ frequentist มักทำงานได้ไม่ดีในกรณีที่ตัวอย่างมีขนาดเล็ก
สถิติแบบเบส์ (Bayesian Statistics)
สถิติแบบเบส์เป็นกระบวนทัศน์ในการวิเคราะห์ข้อมูลอีกลักษณะหนึ่งที่มีความแตกต่างจากไปกระบวนทัศน์แบบดั้งเดิมที่นักวิเคราะห์ข้อมูลทั่วไปคุ้นเคยกัน ซึ่งทำให้สถิติแบบเบส์มีจุดเด่นและสามารถแก้ไขข้อจำกัดที่พบในสถิติแบบ frequentist ได้ ความแตกต่างดังกล่าวเริ่มตั้งแต่ข้อตกลงเบื้องต้นของสถิติแบบเบส์ที่แตกต่างจากสถิติแบบดั้งเดิมอย่างมาก กล่าวคือ
(1) พารามิเตอร์เป็นตัวแปรสุ่ม สถิติแบบเบส์มองว่าพารามิเตอร์เป็นค่าที่ผู้วิเคราะห์สนใจจะศึกษาและต้องการทราบอย่างแท้จริง แต่อย่างไรก็ตามเนื่องจากผู้วิเคราะห์ไม่ทราบข้อมูลประชากร การคาดการณ์หรือความเชื่อ (belief) เกี่ยวกับพารามิเตอร์ดังกล่าวจึงมีความไม่แน่นอน สภาพดังกล่าวจึงสมเหตุสมผลที่จะกำหนด probability statement ให้กับพารามิเตอร์ในโมเดล
จากข้อกำหนดนี้ผลการประมาณค่าพารามิเตอร์ที่ได้จากวิธีการแบบเบส์จึงจะไม่ได้ค่าประมาณแบบค่าเดียว แต่จะได้เป็นการประมาณการแจกแจงของพารามิเตอร์แทน เช่น การประมาณค่าเฉลี่ยผลสัมฤทธิ์ของนักเรียน ผลการวิเคราะห์แบบเบส์จะไม่ได้ให้เป็นค่าประมาณเพียงค่าเดียวคือ 5.84 คะแนน ดัง output ข้างต้น แต่จะให้เป็นการแจกแจงของพารามิเตอร์ค่าเฉลี่ยผลสัมฤทธิ์ทางการเรียน ซึ่งเรียกว่าการแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) และสารสนเทศเกี่ยวกับค่าเฉลี่ยของผลสัมฤทธิ์ดังกล่าวจะถูกสกัดออกมาจากการแจกแจงนี้
(2) ข้อมูลตัวอย่างเป็นค่าคงที่ ถึงแม้ว่าในความเป็นจริงตัวอย่างที่สุ่มมาจากประชากรจะมีความเป็นไปได้ที่หลากหลาย และมีความไม่แน่นอนในการได้มาซึ่งตัวอย่างแต่ละชุด อย่างไรก็ตามเมื่อดำเนินการสุ่มตัวอย่างมาทำการวิเคราะห์แล้ว ความเป็นไปได้จำนวนมากก่อนการสุ่มตัวอย่างจะเหลือเพียงความเป็นไปได้เดียว สถิติแบบเบส์จึงถือว่าข้อมูลตัวอย่างเป็นค่าคงที่ ไม่มีการกำหนด probability statetment ให้กับข้อมูลตัวอย่างนี้
สถิติแบบเบส์สามารถนำมาใช้เพื่อแก้ไขข้อจำกัดที่เกิดขึ้นจากสถิติแบบ frequentist ได้ดังนี้
จากข้อกำหนดเบื้องต้นข้างต้น การอนุมานเชิงสถิติด้วยสถิติแบบเบส์จะดำเนินการผ่านการแจกแจงความน่าจะเป็นของพารามิเตอร์โดยตรงที่เรียกว่า การแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) การแจกแจงดังกล่าวเป็นเครื่องมือสำคัญในการอนุมานเชิงสถิติแบบเบส์ แทน sampling distribution ในการอนุมานเชิงสถิติแบบดั้งเดิม ซึ่งจุดนี้ทำให้สถิติแบบเบส์สร้างความแตกต่างอย่างมากเมื่อเปรียบเทียบกันสถิติแบบดั้งเดิม โดยจะเห็นว่าการประมาณค่าพารามิเตอร์แบบเบส์ไม่ได้อิงกับค่าสถิติเพียงค่าเดียวเหมือนกับสถิติแบบดั้งเดิม แต่อิงกับการแจกแจงความน่าจะเป็นภายหลัง
การที่สถิติแบบเบส์ใช้การแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์เป็นเครื่องมือหลักในการอนุมานเชิงสถิติ จึงทำให้การวิเคราะห์ทำได้อย่างตรงไปตรงมาโดยไม่ต้องอาศัยทฤษฎีหรือวิธีการที่ซับซ้อนเหมือนกับสถิติแบบดั้งเดิม กล่าวคือผู้วิเคราะห์ใช้เพียงสถิติพื้นฐาน เช่น ค่าเฉลี่ย มัธยฐาน ฐานนิยม ส่วนเบี่ยงเบนมาตรฐาน หรือเปอร์เซ็นไทล์ เป็นเครื่องมือในการสรุปสารสนเทศเกี่ยวกับพารามิเตอร์ที่สนใจออกมาจากการแจกแจงความน่าจะเป็นภายหลังที่ประมาณได้ ก็เพียงพอแล้ว นอกจากนี้การแปลผลเพื่ออนุมานเกี่ยวกับพารามิเตอร์ยังทำได้ง่ายและตรงไปตรงมา และได้ข้อสรุปโดยตรงไปยังปริภูมิของพารามิเตอร์ที่เป็นเป้าหมาย แตกต่างจากข้อสรุปที่ได้จากสถิติแบบดั้งเดิมที่เป็นข้อสรุปเกี่ยวกับการแจกแจงของค่าสถิติซึ่งอยู่บนปริภูมิตัวอย่างเท่านั้น
การประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ในโมเดล จะใช้สารสนเทศสองส่วนด้วยกัน ส่วนแรกเป็นสมมุติฐานหรือความเชื่อเบื้องต้น (prior belief) เกี่ยวกับพารามิเตอร์ที่สนใจ ซึ่งกำหนดอยู่ในรูปการแจกแจงความน่าจะเป็นที่เรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์ (prior distribution) และส่วนที่สองคือสารสนเทศที่ได้จากข้อมูลตัวอย่าง ซึ่งแตกต่างจากสถิติแบบดั้งเดิมที่จะใช้สารสนเทศจากข้อมูลตัวอย่างเท่านั้น
ลองพิจารณาตัวอย่างต่อไปนี้
สมมุติว่าต้องการคาดการณ์คะแนนสอบของนายบูม ที่มีความเป็นไปได้ 4 แบบดังรูป
ก่อนการเก็บรวบรวมข้อมูลจริงในชั้นเรียนของนายบูม สมมุติว่าผู้วิเคราะห์มีข้อมูลในอดีตที่เกี่ยวข้องกับคะแนนสอบของนายบูมดังนี้
ไม่มีข้อมูลใด ๆ ที่เป็นประโยชน์ (ให้น้ำหนักกับความเป็นไปได้ A1, A2, A3 และ A4 อย่างไร?)
หากทราบเพิ่มเติมว่าการสอบที่ผ่าน ๆ มาโดยเฉลี่ยนายบูมสอบได้คะแนนคิดเป็นอันดับที่ 4 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 จะเหมือนเดิมมั้ย?)
หากทราบเพิ่มเติมอีกว่า การสอบที่ผ่าน ๆ มา แก้วซึ่งเป็นเพื่อนในห้องเรียนของนายบูม สอบได้คะแนน 38 คะแนน คิดเป็นคะแนนที่น้อยที่สุดในชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)
หากทราบอีกว่า ในการสอบที่ผ่าน ๆ มา นิดสอบได้คะแนนโดยเฉลี่ยคิดเป็น 89 คะแนน ซึ่งเป็นอันดับที่ 3 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)
การแจกแจงก่อนหน้าของคะแนนสอบนายบูม
ผู้วิเคราะห์สามารถปรับสารสนเทศในการแจกแจงความน่าจะเป็นก่อนหน้าในข้างต้นได้ หากมีการเก็บรวบรวมข้อมูลตัวอย่างเพิ่มเติม เช่น
สมมุติว่าการแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมเป็นแบบ uniform และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป
อีกตัวอย่างหนึ่ง สมมุติว่า การแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมมีน้ำหนักสูงที่ A3 และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป
ทฤษฎีบทของเบส์ (Bayes’ theorem)
คำถามที่เกิดขึ้นจากตัวอย่างในข้างต้นคือ จะสามารถจัดสรรหรือปรับปรุงน้ำหนักให้กับแต่ละความเป็นไปได้ในการแจกแจงความน่าจะเป็นก่อนหน้าอย่างไร เมื่อมีข้อมูลเพิ่มเติม คำตอบคือ เราสามารถใช้ทฤษฎีบทของเบส์เพื่อดำเนินการดังกล่าวได้ ทฤษฎีบทของเบส์มีรากฐานมาจากการคำนวณความน่าจะเป็นแบบมีเงื่อนไขดังนี้ \(P(A|B)=P(A \cap B)/P(B)\)
ดังนั้นหากกำหนดให้ \(\theta\) คือพารามิเตอร์ภายในโมเดล และ \(\bf{x}\) คือข้อมูลตัวอย่าง การแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ที่ต้องการคือ \(p(\theta|\bf{x})\) ซึ่งสามารถคำนวณได้จาก
\(p(\theta|\bf{x})=\frac{p(\theta,\bf{x})}{p(\bf{x})}\)
จากสมการในรูปข้างต้น เราสามารถปรับรูปสมการใหม่ได้ดังนี้
\(p(\theta|\bf{x})p(\bf{x})=p(\theta,\bf{x})\) ——- (1)
ในทางกลับกันจากสมการความน่าจะเป็นแบบมีเงื่อนไขจะได้ว่า \(p(\bf{x}|\theta)=p(\theta, \bf{x})/p(\bf{\theta})\) ดังนั้น
\(p(\bf{x}|\theta)p(\bf{\theta})=p(\theta,\bf{x})\) ——- (2)
จากสมการที่ (1) กับ (2) จะได้ว่า
\(p(\theta|\bf{x})p(\bf{x})=p(\bf{x}|\theta)p(\bf{\theta})\)
ดังนั้นจะได้ว่าการแจกแจงความน่าจะเป็นภายหลังสามารถคำนวณได้จาก
เมื่อ \(p(\bf{\theta}|\bf{x})\) คือการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ \(\theta\) จะเห็นว่าเป็นการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ดังกล่าวเมื่อกำหนดข้อมูลตัวอย่าง \(\bf{x}\) ส่วน \(p(\bf{x}|\theta)\) คือฟังก์ชันภาวะความควรจะเป็น \(p(\theta)\) คือฟังก์ชันความน่าจะเป็นก่อนหน้าของพารามิเตอร์ (prior distribution) และ \(p(\bf{x})\) คือค่าคงที่
สมการข้างต้นจึงสามารถเขียนใหม่ได้ดังนี้
จากสมการในข้างต้นจะเห็นว่าการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ เกิดจากการนำสารสนเทศสองส่วนเข้ามาประมวลผลร่วมกัน ได้แก่ สมมุติฐานหรือความเชื่อเบื้องต้น (prior belief) เกี่ยวกับพารามิเตอร์ในโมเดล ที่กำหนดผ่านการแจกแจงความน่าจะเป็นก่อนหน้า และสารสนเทศจากข้อมูลตัวอย่างภายในฟังก์ชันภาวะความควรจะเป็น
ตัวอย่างการอนุมานเชิงสถิติแบบเบส์
โดยปกติขั้นตอนสำคัญของการอนุมานเชิงสถิติแบบเบส์มี 5 ขั้นตอน ได้แก่
การระบุขอบเขตด้านตัวแปร
กำหนดโมเดลของค่าสังเกต
กำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์
ประมาณการแจกแจงความน่าจะเป็นภายหลัง
วิเคราะห์ผลลัพธ์ที่ได้จากการแจกแจงความน่าจะเป็นภายหลัง
หัวข้อนี้จะแสดงตัวอย่างการอนุมานเชิงสถิติแบบเบส์ โดยใช้ตัวอย่างง่าย ๆ ดังนี้
กำหนดตัวแปรในการวิเคราะห์
สมมุติว่าผู้วิจัยต้องการวิเคราะห์ความลำเอียงในการออกหน้าหัวหรือก้อยของเหรียญที่ผลิตขึ้นมารุ่นหนึ่ง การวิเคราะห์นี้ตัวแปรที่สนใจจึงเป็นผลลัพธ์ที่สังเกตได้จากการโยนเหรียญในแต่ละครั้ง กำหนดให้ \(y_i\) แทนค่าสังเกตที่ได้จากการโยนเหรียญในครั้งที่ \(i\) จะได้ว่า
\(y_i=\left\{\begin{matrix} 1 & ; H\\ 0 & ; T \end{matrix}\right.\)
กำหนดโมเดลของค่าสังเกต
เมื่อกำหนดตัวแปรที่สนใจแล้วจะพบว่าตัวแปรที่สนใจมีเพียงตัวแปรตาม 1 ตัว โดยมีค่าสังเกตที่เป็นไปได้เพียง 2 ค่าได้แก่ H หรือ T ซึ่งแทนด้วย 1 กับ 0 ตามลำดับ โมเดลของค่าสังเกต \(y_i\) นี้จึงสามารถกำหนดได้ด้วยการแจกแจงความน่าจะเป็นแบบ Bernoulli กล่าวคือ ความน่าจะเป็นที่จะออกหน้าหัวมีค่าเท่ากับ \(p(y_i=1|\theta)=\theta\) และความน่าจะเป็นที่ออกหน้าก้อยมีค่าเท่ากับ \(p(y_i=0|\theta)=1-\theta\) ซึ่งเขียนเป็นสมการรวมได้ดังนี้
\(p(y_i|\theta)=\theta^{y_i}(1-\theta)^{1-y_i}\) เมื่อ \(i = 1,2,3,...,n\)
จากโมเดลในข้างต้นจะเห็นว่าการเกิดค่าสังเกต \(y_i\) ขึ้นอยู่กับพารามิเตอร์ \(\theta\) ที่มีค่าอยู่บนช่วง \([0,1]\) หากพารามิเตอร์ดังกล่าวมีค่าเท่ากับ \(\theta=0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวและก้อยมีค่าเท่ากัน แต่ถ้าหากพารามิเตอร์ \(\theta>0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวมีมากกว่า จากความหมายดังกล่าวทำให้สามารถมองได้ว่าพารามิเตอร์ \(\theta\) เป็นพารามิเตอร์ความลำเอียงของเหรียญ หากประมาณพารามิเตอร์ดังกล่าวได้จะสามารถสรุปความลำเอียงของเหรียญได้
เนื่องจากเราไม่ได้เก็บรวบรวมข้อมูลเพียงค่าสังเกตเดียว แต่เราทำการทดลองซ้ำ ๆ จำนวน \(n\) ครั้ง เนื่องจากการทดลองแต่ละครั้งเป็นอิสระซึ่งกันและกัน ทำให้โมเดลของชุดค่าสังเกตหรือเรียกอย่างเป็นทางการว่า ฟังก์ชันภาวะความควรจะเป็นดังนี้
\(p(\bf{y}|\theta)=p(y_1|\theta)p(y_2|\theta)...p(y_n|\theta)\)
ซึ่งมีค่าเท่ากับ
\(p(\bf{y}|\theta)=\theta^{\sum_{i=1}^ny_i}(1-\theta)^{n-\sum_{i=1}^ny_i}\)
ในฟังก์ชันภาวะความควรจะเป็นนี้จะให้ค่าความน่าจะเป็นหรือค่าความหนาแน่นของข้อมูลตัวอย่าง \(\bf{y}\) บนแต่ละค่าที่เป็นไปได้ของพารามิเตอร์ \(\theta\) เนื่องจากข้อมูลตัวอย่างเป็นค่าคงที่ ดังนั้นค่าความหนาแน่นนี้จึงเปลี่ยนแปลงไปตามค่าของพารามิเตอร์ \(\theta\) โดยหากค่าความหนาแน่นนี้มีค่าสูง ณ ค่า \(\theta\) ค่าใด นั่นหมายถึงว่าพารามิเตอร์ค่านั้นทำให้โมเดลของค่าสังเกตกับค่าสังเกตจริงมีความสอดคล้องกันสูง ข้อสังเกตนี้ถูกนำไปใช้เพื่อหาค่าประมาณด้วยวิธีภาวะความควรจะเป็นสูงสุด (maximum likelihood: ML) แต่ในวิธีการแบบเบส์จะนำสารสนเทศส่วนนี้ไปประมวลร่วมกับความเป็นไปได้ของพารามิเตอร์ในการแจกแจงความน่าจะเป็นเบื้องต้นก่อน
กำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์
การกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) สามารถทำได้หลายลักษณะ ทั้งแบบไม่ต่อเนื่อง และแบบต่อเนื่อง ขึ้นอยู่กับขอบเขตของการวิเคราะห์ ข้อมูลในอดีตที่ใช้สนับสนุน และธรรมชาติของพารามิเตอร์นั้น ในตัวอย่างข้างต้นพารามิเตอร์ \(\theta\) ต้องมีค่าอยู่บนช่วง \([0,1]\)
เพื่อให้ตัวอย่างนี้ทำความเข้าใจได้ง่าย จึงเลือกกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) เป็นการแจกแจงแบบไม่ต่อเนื่อง โดยมีค่าที่เป็นไปได้เท่ากับ 0.0, 0.1, 0.2, 0.3, …., 1.0 และมีการแจกแจงของความน่าจะเป็นก่อนหน้าดังรูป
# prior distribution
theta<-seq(0,1,0.1)
p.theta<-c(0.01,0.04,0.075,0.1,0.15,0.25,0.15,0.1,0.075,0.04,0.01)
plot(theta,p.theta,
xlab=expression(theta),
ylab=expression(p(theta)),
type="h", lwd=2, col="orange",
main="Prior Distribution")
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")
การประมาณการแจกแจงความน่าจะเป็นภายหลัง
จากการทดลองโยนเหรียญจำนวน 10 ครั้งพบว่า ออกหน้าหัวจำนวน 6 ครั้ง ดังนั้นจะได้ว่า
\(p(\bf{y}|\theta)=\theta^{6}(1-\theta)^4\)
จากขอบเขตของ \(theta\) ที่กำหนดในการแจกแจงความน่าจะเป็นเบื้องต้นจะได้ว่าฟังก์ชันภาวะความควรจะเป็น มีลักษณะดังรูป
# likelihood function
L<-theta^6*(1-theta)^4
plot(theta,L,
xlab=expression(theta),
ylab="p(y|theta)",
type="h", lwd=2, col="orange",
main="Likelihood Function")
points(theta,L, type="p",pch=16, cex=2, col="orange")
จากทฤษฎีบทของเบส์เราสามารถรวมสารสนเทศจากการแจกแจงความน่าจะเป็นเบื้องต้น กับฟังก์ชันภาวะความควรจะเป็นได้ดังนี้
# marginal y
p.y<-sum(p.theta*L)
par(mfrow=c(1,2))
# posterior of theta given y
post.theta<-L*p.theta/p.y
plot(theta,post.theta,
xlab=expression(theta),
ylab="p(theta|y)",
type="h", lwd=2, col="blue",
ylim=c(0,0.5),
main="Posterior Distribution")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")
# prior and posterior in the same plot
plot(theta,p.theta,
xlab=expression(theta),
ylab="Probability",
type="l", lwd=2, col="orange",
ylim=c(0,0.5))
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")
points(theta,post.theta,
type="l", lwd=2, col="blue")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")
legend(-0.05,0.48, legend=c("Prior","Posterior"), col=c("orange","blue"), lty=1, bty="n", lwd=2)
การอนุมานเชิงสถิติจากการแจกแจงความน่าจะเป็นภายหลัง
เมื่อได้การแจกแจงความน่าจะเป็นภายหลังที่ต้องการแล้ว ผู้วิเคราะห์สามารถวิเคราะห์การแจกแจงดังกล่าวเพื่ออนุมานเกี่ยวกับพารามิเตอร์ความลำเอียงที่ต้องการศึกษาได้ โดยสามารถทำได้หลายลักษณะ ทั้งการประมาณค่าแบบจุด การประมาณค่าแบบช่วง และการทดสอบสมมุติฐาน รายละเอียดดังนี้
(1) การประมาณค่าแบบจุด
การประมาณค่าแบบจุดสามารถทำได้โดยใช้สถิติพื้นฐานสรุปแนวโน้มสู่ส่วนกลางของค่าพารามิเตอร์ จากการแจกแจงความน่าจะเป็นภายหลัง เช่น ค่าเฉลี่ย มัธยฐาน และฐานนิยม นอกจากนี้ยังสามารถคำนวณค่าส่วนเบี่ยงเบนมาตรฐานเพื่อประเมินระดับความน่าเชื่อถือของค่าประมาณแบบจุดดังกล่าวได้อีกด้วย จากการแจกแจงความน่าจะเป็นภายหลังข้างต้น จะได้ว่า
# mean = expected value of theta
mean.theta<-sum(post.theta*theta)
sd.theta<-sqrt(sum(post.theta*(theta-mean.theta)^2))
mean.theta
## [1] 0.5540441
sd.theta
## [1] 0.1145753
(2) การประมาณค่าแบบช่วง
ช่วงการประมาณที่ได้จากวิธีแบบเบส์จะเรียกว่า ช่วงความน่าเชื่อถือ (credible interval) การประมาณค่าพารามิเตอร์แบบช่วงจากการแจกแจงความน่าจะเป็นภายหลังสามารถทำได้หลายวิธีการ วิธีการง่าย ๆ คือการใช้ค่าเฉลี่ยบวกลบด้วยส่วนเบี่ยงเบนมาตรฐานของพารามิเตอร์ อีกวิธีการหนึ่งที่มีประสิทธิภาพมากกว่าคือการหาช่วงแบบ highest density interval (HDI) กล่าวคือเป็นช่วงการประมาณที่มีความหนาแน่นมากที่สุดที่ทำให้ค่าความน่าจะเป็นของพารามิเตอร์ภายในช่วงดังกล่าวเท่ากับค่าที่กำหนดเช่น ช่วง 95% HDI ของพารามิเตอร์ \(\theta\) ในข้างต้นจะมีค่าประมาณ \([0.4,0.8]\)
ช่วงดังกล่าวมีความแตกต่างจาก 95% confidence interval อย่างไร??
(3) การทดสอบสมมุติฐาน
การทดสอบสมมุติฐานด้วยวิธีการแบบเบส์มีความยืดหยุ่นสูงมาก และสามารถทำได้หลายวิธีการ วิธีการแรกเรียกว่า maximum a posteriori (MAP) test มีรายละเอียดดังนี้
กำหนดให้ \(H_0\) และ \(H_1\) เป็นสมมุติฐานที่ต้องการเปรียบเทียบ ในกรณีนี้จะเห็นว่ามีความไม่แน่นอนว่าสมมุติฐานใดเป็นสมมุติฐานที่ถูกต้องกันแน่ ดังนั้นจึงสามารถกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าของสมมุติฐานทั้งสองได้ ดังนี้
\(p(H_0) = p_0\) และ \(p(H_1)=p_1\) โดยที่ \(p_0+p_1=1\)
และกำหนดให้ \(p(\bf{y}|H_0)\) กับ \(p(\bf{y}|H_1)\) คือฟังก์ชันภาวะความควรจะเป็นบนสมมุติฐานทั้งสอง ซึ่งจากทฤษฎีบทของเบส์จะได้ว่า
\(p(H_0|\bf{y})=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y})}\)
\(p(H_1|\bf{y})=\frac{p(\bf{y}|H_1)p(H_1)}{p(\bf{y})}\)
การตัดสินใจว่าควรเลือกเชื่อสมมุติฐานใด สามารถทำได้โดยดูจากอัตราส่วนที่เรียกว่า posterior odd ซึ่งมีค่าเท่ากับ
\(\frac{p(H_0|\bf{y})}{p(H_1|\bf{y})}=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y}|H_1)p(H_1)}\)
โดยหากอัตราส่วนข้างต้นมีค่ามากกว่า 1 แสดงว่า สมมุติฐาน \(H_0\) มีแนวโน้มน่าเชื่อถือมากกว่า ในทางกลับกันหากอัตราส่วนดังกล่าวมีค่าน้อยกว่า 1 แสดงว่าสมมุติฐาน \(H_1\) มีแนวโน้มน่าเชื่อถือมากกว่า
จากตัวอย่างโยนเหรียญ หากต้องการทดสอบว่าเหรียญมีความเที่ยงตรงหรือไม่ อาจกำหนดสมมุติฐานเป็นดังนี้
\(H_0: \theta = 0.5\) vs \(H_1: \theta>0.5\)
prior.H<-c(0.5,0.5) # prior for H0 and H1 respectively.
# calculate likelihood
L.H0<-0.5^6*(1-0.5)^4
theta.H1<-theta[theta>0.5]
L.H1<-sum(theta.H1^6*(1-theta.H1)^4)
# calculate posterior odd
(prior.H[1]*L.H0)/(prior.H[2]*L.H1)
## [1] 0.3727444
อีกวิธีการหนึ่งที่สามารถทำได้เรียกว่า minimum cost hypothesis test รายละเอียดมีดังนี้
กำหนดให้
- \(C_{10}\) คือ cost ของการเลือก \(H1\) โดยที่ \(H_0\) ถูกต้อง (cost ของการเกิด type I error)
- \(C_{01}\) คือ cost ของการเลือก \(H_0\) โดยที่ \(H_1\) ถูกต้อง (cost ของการเกิด type II error)
เกณฑ์การพิจารณาจะใช้ค่าความเสี่ยงภายหลัง (posterior risk) ดังนี้
\(\frac{p(H_0|\bf{y})C_{10}}{p(H_1|\bf{y})C_{01}}\)
จากตัวอย่างโยนเหรียญสมมุติว่า \(C_{10}=5C_{01}\)
c10<-5
c01<-1
# posterior risk of accepting H0
(prior.H[2]*L.H1)*c01
## [1] 0.001309962
# posterior risk of accepting H1
(prior.H[1]*L.H0)*c10
## [1] 0.002441406
# risk ratio
(prior.H[1]*L.H0)*c10/(prior.H[2]*L.H1)*c01
## [1] 1.863722
Monte Carlo Markov Chain (MCMC)
ตัวอย่างที่แสดงให้ดูในหัวข้อข้างต้นเป็นตัวอย่างที่ง่ายมากทำให้การคำนวณการแจกแจงความน่าจะเป็นภายหลังสามารถคำนวณมือได้โดยสะดวก นอกจากนี้ยังสามารถพิสูจน์ในรูปของสูตรปิดทางคณิตศาสตร์ได้ด้วย อย่างไรก็ตามในสถานการณ์ทั่วไปการคำนวณการแจกแจงความน่าจะเป็นภายหลังดังกล่าวมักมีความซับซ้อนและไม่สามารถคำนวณมือได้ โดยเฉพาะในโมเดลพหุระดับ ปัจจุบันวิธีการที่ถูกนำมาแทนการคำนวณทางคณิตศาสตร์คือการใช้ลูกโซ่มาร์คอฟมอนติคาร์โล (MCMC) วิธีการนี้ถูกออกแบบให้สุ่มตัวอย่างพารามิเตอร์ภายใต้กระบวนการลูกโซ่ของมาร์คอฟ ซึ่งรับประกันว่าเมื่อจำนวนรอบของการสุ่มตัวอย่างมีจำนวนมากเพียงพอ ตัวอย่างของพารามิเตอร์ดังกล่าวจะถูกสุ่มมาจากการแจกแจงของประชากรที่เป็นการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ตามต้องการ ในการทำงานผู้วิเคราะห์จะใช้ตัวอย่างสุ่มที่ถูกตรวจสอบคุณสมบัติแล้ว มาประมาณการแจกแจงความน่าจะเป็นภายหลัง การอนุมานเชิงสถิติต่าง ๆ ที่ได้กล่าวไว้บ้างก่อนหน้านี้สามารถทำได้โดยตรงผ่านการวิเคราะห์ตัวอย่างของพารามิเตอร์ดังกล่าว
ปัจจุบันโปรแกรมสำเร็จรูปที่ช่วยทำ MCMC ให้กับผู้วิเคราะห์มีหลายตัว ได้แก่ BUGS, JAGS, Stan รวมทั้ง R อีกหลาย package นอกจากนี้โปรแกรม Mplus และ MLWins ยังมีโมดูลสำเร็จรูปสำหรับประมาณค่าพารามิเตอร์แบบเบส์สำหรับโมเดลพหุระดับอีกด้วย
บทเรียนนี้จะทำ MCMC ด้วย 2 วิธีการ วิธีการแรกคือการทำ MCMC บนโปรแกรม JAGS ที่สั่งการทำงานบนโปรแกรม R และวิธีการที่สองคือการทำ MCMC ด้วย package-brms ซึ่งเป็น high-level API ของโปรแกรม Stan โดยถูกพัฒนาขึ้นมาสำหรับวิเคราะห์โมเดลพหุระดับโดยเฉพาะ จุดเด่นของ package-brms คือใช้หลักภาษาที่ง่ายเกือบจะเหมือนกับ package-lme4 ที่ได้เรียนไปในบทเรียนก่อนหน้า นอกจากนี้ยังสามารถวิเคราะห์ได้หลายโมเดลตั้งแต่โมเดลแบบตัวแปรตามตัวเดียวไปถึงตัวแปรตามหลายตัว และมีขอบเขตครอบคลุมไปถึง generalized linear model แบบพหุระดับอีกด้วย จุดเด่นอีกประการหนึ่งของ package-brms คืออิงกับโปรแกรม Stan ที่ใช้อัลกอริทึม Hamiltonian Monte Carlo ร่วมกับอัลกอริึม No-U-Turn Sampler (NUTS) ซึ่งมีประสิทธิภาพสูงกว่า Gibb sampler ที่ใช้ใน BUGS และ JAGS อัลกอริทึมดังกล่าวจะลู่เข้าสู่สถานะคงที่ (steady state) ซึ่งก็คือได้ตัวอย่างที่สุ่มจากการแจกแจงความน่าจะเป็นภายหลังไวจึงเหมาะกับโมเดลซับซ้อนที่มีพารามิเตอร์จำนวนมาก นอกจากนี้อัลกอริทึมดังกล่าวยังไม่ต้องการการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าแบบวงศ์คู่สังยุค (conjugacy prior) เหมือนกับอัลกอริทึม Gibb sampler อีกด้วย จึงทำให้มีความยืดหยุ่นในการวิเคราะห์ที่สูงขึ้นมาก อย่างไรก็ตามข้อจำกัดของ Stan คือการสุ่มตัวอย่างในแต่ละรอบจะใช้เวลาที่มากกว่า BUGS และ JAGS (Bürkner, 2017)
JAGS (Just Another Gibb Sampler)
JAGS สามารถรันได้บน 3 platforms ได้แก่ Mac, Windows และ Linux บทเรียนนี้จะควบคุมโปรแกรม JAGS ผ่านโปรแกรม R โดยใช้ package-rjags ดังนั้นผู้วิเคราะห์จะใช้โปรแกรม R ในการทำงานส่วนใหญ่ ท้ังการเตรียมข้อมูล และการวิเคราะห์ผลลัพธ์จากลูกโซ่มาร์คอฟ ส่วนที่จะเขียนเป็นภาษาของ JAGS คือส่วนการระบุโมเดลการวิเคราะห์เท่านั้น การกำหนดโมเดลใน JAGS ใข้หลักภาษาเดียวกับ BUGS โดยทั่วไปโมเดลการวิเคราะห์แบบเบส์ในโปรแกรม JAGS จะมีส่วนประกอบ 3 ส่วน ได้แก่
likelihood function(s)
prior distribution
deterministic model
prior และ likelihood function เป็นส่วนประกอบที่เรียกว่า stochastic node ซึ่งนิยามโดยใช้สัญลักษณ์ ~ ส่วน deterministic model เป็นส่วนค่าคงที่หรือส่วนโมเดลเชิงคณิตศาสตร์ที่ใช้แสดงความสัมพันธ์ระหว่างตัวแปรการนิยามส่วนนี้จะใช้สัญลักษณ์ <-
ตัวอย่างด้านล่างแสดงการเขียนโมเดลการวิเคราะห์การถดถอยอย่างง่ายแบบเบส์ เพื่อวิเคราะห์ความสัมพันธ์ระหว่างผลสัมฤทธิ์ทางการเรียน (ACH) กับความตั้งใจเรียนของนักเรียน (ATT)
กำหนดให้ \(y\) แทน ACH และ \(x\) แทน ATT เนื่องจากตัวแปรตามในโมเดลเป็นตัวแปรต่อเนื่อง ในกรณีนี้จึงกำหนดให้
\(y_i \sim N(\mu_i, \sigma^2)\)
สังเกตว่าการแจกแจงของค่าสังเกตในข้างต้นมีค่าเฉลี่ยที่ห้อย \(i\) ทั้งนี้เป็นเพราะเรากำลังวิเคราะห์การถดถอย ดังนั้นค่าเฉลี่ยของ \(y\) จึงมีการเปลี่ยนแปลงไปตามค่าของตัวแปรอิสระ \(x_i\) ค่าเฉลี่ยดังกล่าวจึงเป็นค่าเฉลี่ยที่มีเงื่อนไข และสามารถเขียนเป็นสมการแสดงความสัมพันธ์ได้ดังนี้
\(\mu_i = \beta_0 + \beta_1x_i\)
โมเดลข้างต้นมีพารามิเตอร์จำนวน 3 ตัว จำแนกเป็นสองกลุ่ม กลุ่มแรกคือพารามิเตอร์อิทธิพลคงที่ (fixed-effect parameter) ได้แก่ \(\beta_0\) และ \(\beta_1\) และกลุ่มที่สองคือพารามิเตอร์ความแปรปรวนของอิทธิพลสุ่ม ได้แก่ \(\sigma^2\)
การแจกแจงความน่าจะเป็นก่อนหน้า
การแจกแจงความน่าจะเป็นก่อนหน้าใช้สะท้อนความเชื่อเบื้องต้นในพารามิเตอร์ต่าง ๆ ของโมเดล ความเชื่อดังกล่าวมาได้จากทั้งผลการศึกษาในอดีต หลักเหตุผล และประสบการณ์ของนักวิจัย นอกจากนี้ยังควรต้องพิจารณาธรรมชาติหรือค่าที่เป็นไปได้ของพารามิเตอร์ประกอบด้วย
โดยทั่วไปการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าสามารถทำได้ 2 ลักษณะ ลักษณะแรกเรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ (non-informative prior distribution) และการแจกแจงความน่าจะเป็นก่อนหน้าแบบให้สารสนเทศ (informative prior distribution)
การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์อิทธิพลคงที่สามารถเลือกใช้ได้หลายการแจกแจง เช่น การแจกแจงแบบ uniform ที่เป็น noninformative prior การแจกแจงแบบปกติ ที่เป็นไปได้ทั้ง noninformative และ informative prior ขึ้นอยู่กับการกำหนดพารามิเตอร์ค่าเฉลี่ยและส่วนเบี่ยงเบนมาตรฐานของการแจกแจง
จากตัวอย่างข้างต้นกำหนดให้
\(\beta_0 \sim N(0,10000)\)
\(\beta_1 \sim N(0,10000)\)
ซึ่งเป็นการแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ
การแจกแจงความน่าจะเป็นก่อนหน้าของส่วนเบี่ยงเบนมาตรฐานหรือความแปรปรวน มีความจำเป็นต้องเลือกการแจกแจงที่ domain ของการแจกแจงมีค่าไม่ติดลบ โดยปกติมักกำหนดให้ความแปรปรวนดังกล่าวมีการแจกแจงแบบแกมมาผกผัน (inverse gamma distribution)
จากตัวอย่างข้างต้น กำหนดให้
\(sigma^2 \sim IG(a,b)\)
โดยที่ \(IG\) คือการแจกแจงแบบ inverse-gamma
Note: ในเชิงทฤษฎีหากกำหนดให้ \(sigma^2 \sim IG(a,b)\) จะได้ว่า
\(p(\sigma^2|a,b) \propto (\sigma^2)^{-a-1}exp(\frac{-b}{\sigma^2})\)
หากกำหนดให้ \(a \rightarrow 0\) และ \(b \rightarrow 0\) แล้วฟังก์ชันความหนาแน่นในข้างต้นจะเขียนใหม่ได้เป็น
\(p(\sigma^2|a,b) \propto \frac{1}{\sigma^2}\) เรียกว่า Jeffery prior ซึ่งเป็น noninformative prior ที่ใช้ได้ตัวหนึ่งของพารามิเตอร์ความแปรปรวน
ดังนั้นในกรณีนี้จึงกำหนดให้
\(sigma^2 \sim IG(0.0001,0.0001)\)
รูปด้านล่างแสดงรายการของการแจกแจงความน่าจะเป็นที่สามารถกำหนดได้ใน BUGS และ JAGS
จากการระบุโมเดลข้างต้นทำให้สามารถเขียน syntax เพื่อระบุโมเดลใน JAGS ได้ดังนี้
model{
for(i in 1:n)
{
# likelihood function
y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
# deteministic model
mu[i]<-b0+b1x[i]
}
# prior distribution
b0~dnorm(0,10^-4)
b1~dnorm(0,10^-4)
tau~gamma(0.0001,0.0001)
sigma2<-1/tau # deteministic model
}
ขั้นตอนการประมวลผลด้วยโปรแกรม JAGS
การประมวลผลด้วย JAGS มี 4 ขั้นตอนหลักได้แก่
เตรียมข้อมูล
ระบุโมเดล
ขั้นตอนนี้สามารถทำได้ดังตัวอย่างข้างต้น และเมื่อระบุโมเดลในโปรแกรม R เรียบร้อยแล้ว ผู้วิเคราะห์จำเป็นต้องเขียนโมเดลดังกล่าวในรูปของไฟล์ .txt โดยใช้คำสั่ง writeLines(modelString, con="model.txt") ดังตัวอย่างด้านล่าง
modelString<-"model{
for(i in 1:n)
{
# likelihood function
y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
# deteministic model
mu[i]<-b0+b1*x[i]
}
# prior distribution
b0~dnorm(0,10^-4)
b1~dnorm(0,10^-4)
tau~dgamma(0.0001,0.0001)
sigma2<-1/tau # deteministic model
}" #Wrap model syntax into String object
writeLines(modelString, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/simReg.txt") # write modelString to model.txt
ส่งข้อมูลและโมเดลไปประมวลผลบน JAGS
ตรวจสอบการลู่เข้าของ MCMC
วิเคราะห์และแปลผล
สมมุติว่าข้อมูลที่ใช้ในการวิเคราะห์เป็นดังรูป โดยเก็บไว้ในตัวแปร dat.reg
head(dat.reg)
## ATT ACH
## 1 30.69309 56.96164
## 2 11.95094 28.41848
## 3 32.41264 42.84734
## 4 20.02710 51.75678
## 5 19.52567 32.15973
## 6 37.61670 57.54749
par(mar=c(4,4,1,1))
plot(ATT,ACH, pch=16, xlab="ATT", ylab="ACH")
การประมวลผลโมเดลข้างต้นผ่านโปรแกรม JAGS จะต้องใช้คำสั่ง 2 ตัวได้แก่ฟังก์ชัน jag.model() และ coda.samples() ที่มีอากิวเมนท์สำคัญดังนี้
setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")
library(rjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)
model<-jags.model(file = "simReg.txt",
data = dataList,
inits = initsList,
n.chains = 3)
update(model, n.iter=2000) # n.burnin
samples<-coda.samples(model,
variable.names=c("b0","b1","sigma2"),
n.iter = 10000,
thin=1)
head(samples)
plot(samples)
summary(samples)
อีก package หนึ่งที่สามารถเรียก JAGS จาก R ได้คือ package-runjags ซึ่งมีจุดเด่นคือสามารถสั่งประมวลแบบคู่ขนาน (parallel) แยกตาม core ของ CPU ได้ในกรณีที่ผู้วิเคราะห์กำหนดให้สุ่มตัวอย่างจากลูกโซ่มาร์คอฟที่มีจำนวนมากกว่า 1 ลูกโซ่ ตัวอย่างคำสั่งเป็นดังนี้
setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")
library(runjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)
runjag.output<-run.jags(method="parallel",
model="simReg.txt",
monitor=c("b0","b1","sigma2"),
data=dataList,
inits=initsList,
n.chains=3,
burnin=2000,
sample=10000, #n.iter
thin=1,
summarise=FALSE,
plots=FALSE)
samples<-as.mcmc.list(runjag.output)
head(samples)
plot(samples)
summary(samples)
การวิเคราะห์ผลลัพธ์ที่ได้จาก MCMC
การวิเคราะห์ผลลัพธ์ที่ได้จาก MCMC มีวัตถุประสงค์ 2 ประการ ประการแรกคือการวิเคราะห์เพื่อตรวจสอบ/วินิจฉัยว่าตัวอย่างที่ได้จาก MCMC ในข้างต้น เป็นตัวอย่างที่สุ่มได้จากการแจกแจงความน่าจะเป็นภายหลังที่เป็นเป้าหมายแล้วหรือไม่ (ลูกโซ่ที่สร้างขึ้นลู่เข้าสู่ posterior distribution เป้าหมายแล้วหรือไม่) การตรวจสอบดังกล่าว
ผลลัพธ์ที่ได้จากอัลกอริทึม MCMC จะเป็นตัวอย่างของพารามิเตอร์ในโมเดลที่มีจำนวนเท่ากับ จำนวนรอบทวนซ้ำ (iteration) ตามที่กำหนด หากจำนวนรอบที่ทวนซ้ำมีมากเพียงพอตัวอย่างดังกล่าวจะเป็นค่าพารามิเตอร์ที่ถูกสุ่มมาจากการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ภายในโมเดล ผู้วิเคราะห์สามารถใข้ข้อมูลตัวอย่างดังกล่าวในการอนุมานเชิงสถิติผ่านการแจกแจงของตัวอย่างพารามิเตอร์ อย่างไรก็ตามก่อนการใช้งานตัวอย่างดังกล่าวผู้วิเคราะห์จำเป็นต้องตรวจสอบคุณสมบัติของตัวอย่างพารามิเตอร์ดังกล่าวก่อน เพื่อรับประกันความถูกต้องของผลการวิเคราะห์ โดยมีประเด็นการตรวจสอบดังนี้
- อัตสหสัมพันธ์ระหว่างตัวอย่างภายในแต่ละลู่โซ่ (autocorrelation) เนื่องจากตัวอย่างของพารามิเตอร์ถูกจำลองขึ้นโดยสุ่มจากกระบวนการมาร์คอฟ ดังนั้นจึงมีความเป็นไปได้ที่ตัวอย่างในแต่ละรอบจะมีความสัมพันธ์ซึ่งกันและกัน หากความสัมพันธ์ดังกล่าวมีค่าสูงมากเกินไปจะทำให้ผลการวิเคราะห์จากลูกโซ่ดังกล่าวมีความคลาดเคลื่อนได้ การแก้ปัญหาดังกล่าวสามารถทำได้ด้วยการลบตัวอย่าง
และผู้วิเคราะห์สามารถใช้ตัวอย่างของพารามิเตอร์ดังกล่าวแทนค่าประมาณของการแจกแจงความน่าจะเป็นภายหลัง
ตัวอย่างโยนเหรียญ (รอบที่ 2)
- จากตัวอย่างที่ 1 จะได้ว่าโมเดลของค่าสังเกตเป็นฟังก์ชันความน่าจะเป็นแบบ Bernuolli ที่มีพารามิเตอร์เป็น \(\theta\) เขียนแทนด้วย
\(y_i \sim Ber(\theta)\) โดยที่ \(i=1,2,...,n\)
- อย่างไรก็ตามในกรณีนี้ผู้วิเคราะห์ต้องการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าแบบต่อเนื่องให้กับพารามิเตอร์ความลำเอียง \(\theta\) เนื่องจากพารามิเตอร์ดังกล่าวต้องมีค่าอยู่บนช่วง \([0,1]\) การแจกแจงความน่าจะเป็นที่เหมาะสมจึงเป็นการแจกแจงแบบบีต้า (beta distribution) ที่มีพารามิเตอร์ตำแหน่ง \(a\) และพารามิเตอร์สเกล \(b\) ดังนั้น
\(\theta \sim Beta(a,b)\)
พารามิเตอร์ทั้งสองตัวของการแจกแจงแบบบีต้าทำให้รูปทรงการแจกแจงมีความแตกต่างกัน ซึ่งสามารถใช้สะท้อนความเชื่อเกี่ยวกับความลำเอียงของเหรียญที่แตกต่างกันได้เป็นอย่างดี รูปด้านล่างแสดงตัวอย่างโค้งความหนาแน่นของการแจกแจงแบบบีต้า เมื่อกำหนดพารามิเตอร์ต่าง ๆ
Kruschke (2012)
จากการกำหนดในข้างต้นจะสามารถเขียน syntax ของโมเดลในโปรแรกม JAGS ได้ดังนี้
model{
# likelihood part
for (i in 1:n)
{
y[i]~dbern(theta)
}
# prior part
theta~beta(a,b)
# deterministic part
a<-3
b<-2
}
การประมวลผลโมเดลข้างต้นด้วยโปรแกรม JAGS บน R จำเป็นต้องมีการ wrap up คำสั่งข้างต้นให้เป็น text file เพื่อส่งออกไปยังโปรแกรม R โดยเขียนคำสั่งดังนี้
coin_model="
model{
# likelihood part
for (i in 1:n)
{
y[i]~dbern(theta)
}
# prior part
theta~dbeta(a,b)
# deterministic part
a<-3
b<-2
}
"
writeLines(coin_model, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt")
ขั้นตอนถัดมาคือการประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ \(\theta\) โดยใช้เทคนิค MCMC ผ่านโปรแกรม JAGS ซึ่งสามารถสั่งงานได้ผ่านฟังก์ชัน jag.model() โดยฟังก์ชันนี้มีอาร์กิวเมนท์ที่สำคัญที่เกี่ยวข้องกับการระบุคุณลักษณะของลูกโซ่มาร์คอฟที่ต้องการสร้างขึ้น ได้แก่
- ข้อมูลค่าสังเกต ในที่นี้สมมุติว่าทำการทดลอง 20 ครั้งได้ผลเป็นดังนี้
y<-rbinom(20,1,0.65)
data.list<-list(y=y, n=20)
n.chains
library(rjags)
fit<-jags.model(file="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt", data=data.list,
n.chains = 3)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 20
## Unobserved stochastic nodes: 1
## Total graph size: 24
##
## Initializing model
samples<-coda.samples(fit,variable.names=c("theta"),n.iter=5000)
plot(samples)
summary(samples)
##
## Iterations = 1:5000
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 0.7191290 0.0878639 0.0007174 0.0007152
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## 0.5340 0.6610 0.7243 0.7830 0.8750